Introduction

A good quality and the right amount of surface water have been crucial for human well-being and life in general ever since (Water Framework Directive 2000). These two aspects: quality and quantity, are determined to a significant degree by the ecosystems surrounding the water bodies (Dingman 2015; Erős et al. 2019; Şandric et al. 2019). The quality of water can be impaired by processes within the catchment area. Pressures exerted by a diversity of human activities (e.g. application of fertilizers, pesticides, removal of vegetation, ploughing, soil surface sealing) are transmitted, and their effects can potentially accumulate in the water bodies (Chapra 1997). On the other hand, vegetation cover can modify these processes, adding to or alleviating the pressures (Vought et al. 1994; Olde Venterink et al. 2006; Conte et al. 2011). The actual ratio of different land use–land cover (LULC) types nowadays needs active decision-making, especially to encounter adverse effects of anthropogenic pressures.

Hydrologic ecosystem services (HES) are mostly regulating services (classification in the Common International Classification of Ecosystem Services (CICES): Regulation and Maintenance), regulating the quantity or the quality of water (Haines-Young and Potschin 2017). In contrast to aquatic ecosystem services (Martin-Ortega et al. 2015), which encompass ecosystem services (ES) related to the aquatic environment that take place within the water body (Brauman 2016), HES are determined by hydrologic processes within the catchment area. Many results of the processes can be seen in the aquatic environment, e.g. in fauna and flora of the water bodies, which can be detected by water quality monitoring, for example, within the Water Framework Directive (WFD). By regulating the movement of water, not just surface and subsurface flow of water is slowed down, thus regulating the accumulation of precipitation (flood regulation or flood control ES), but also influencing other transportation processes (control of soil erosion, filtration of pollutants, moderating their spread/distribution towards downstream water).

While the implementation of the ES concept suggests a complex, objective valuation of natural capital throughout the landscape (Millennium Ecosystem Assessment 2005), it is not at all evident how the different aspects/services addressed in an ES assessment should be reconciled and integrated (Sutherland et al. 2018). Mapping, modelling and quantifying ES on watersheds can significantly support decision-making in land use conflict areas based on multi-faceted analyses (Crossman et al. 2013; Nedkov et al. 2015; Burkhard and Maes 2017). For management and planning, spatially explicit decisions need to be made; thus, the mapping of different ES in an aggregated way is of special importance.

Changes in land use are mostly incentivized by the subvention system or Common Agricultural Policy (CAP) (Tokarick 2005; Hemming et al. 2018), while for specific mitigating and restoration measures, funding is often very limited. Thus, for a successful management not just good decision-making practice is needed, but the localization of critical sites for most effective intervention is also essential. Trade-offs between different ES need to be taken into account in order to arrive at the best possible solutions (Raudsepp-Hearne et al. 2010).

River Zala is the largest tributary stream to Lake Balaton, which is of outstanding economic and social value in Hungary, especially as a national and international holiday resort. Thus, the quantity and quality of its water are of great importance, in order to ensure good ecological status of the lake as well as uninhibited recreation (Pomucz and Csete 2020).

The goal of our research was (i) to model and map selected regulating HES in a major catchment area of Lake Balaton, (ii) to suggest a simple procedure for aggregating results from ES mapping, (iii) to localize critical and outstanding sites based on the aggregated results, (iv) to design minimal changes in a future land use scenario based on the aggregated ES results, focusing on critical sites and (v) to calculate the differences in the landscape’s capacity to provide HES under the hypothesized future land use.

Materials and methods

Study area

We dealt with the Zala River catchment area. Zala River is the largest tributary of Lake Balaton, the largest shallow lake in Central Europe. Thus, Zala River’s mouth was modified in the 1980s and the Kis-Balaton Water Protection System (KBWPS) was created. The aim was to improve the deteriorating water quality of Lake Balaton by artificially restoring the buffering function of former wetlands, which retain sediment and nutrients from the Zala River (Pomogyi 1993; Tátrai et al. 2000). As the constructed wetland system is dominated by different water management and complex, dynamic water quality processes, we only considered the catchment area upstream of the KBWPS. The studied area is almost 1530 km2 (Fig. 1).

Fig. 1
figure 1

Study area

Full size image

The LULC conditions of the study area were determined from the National Ecosystem Map, developed within the Mapping and Assessment of Ecosystems and their Services in Hungary (MAES-HU), as part of the implementation of the EU Biodiversity Strategy to 2020 in Hungary (European Comission 2011; Kovács-Hostyánszki et al. 2019; Tanács et al. 2019). Considering first-level categories (the map was implemented with three-level structure), LULC is dominated by forests (42%), while croplands (vineyards, orchards and arable lands) also have a high coverage (38%). The urban areas (residential, industrial and urban areas) (11%), grasslands (6%) and water bodies (including rivers, lakes and wetlands) (3%) make up the remaining land. The most populous town in the study area is Zalaegerszeg (with about 60,000 inhabitants).

The catchment is dominated by luvisol soils. According to the latest soil maps (Pásztor et al. 2017), the texture in the upper 30 cm layer is covered almost half of the area with loam (50%) and the other half of the basin is topped by silty loam (34%) and sandy loam (16%).

The climate is moderately cool and moderately humid. The average annual cloud cover is 55–65%, with the degree of cloud cover decreasing from west to south. The number of hours with sunshine is 1900–2000 h per year. The spatial distribution of long-term average of precipitation time series shows an increasing trend from the east to the west (annually 660 mm in the north-east corner and 800 mm in the western corner). The mean air temperature follows a reverse trend, being typically 1.0–1.5 °C cooler in the west than in the north-east (based on the Second Hungarian River Basin Management Plan).

The riverbed of Zala River is around 7 to 20 m wide, and the water depth is between 0.5 and 2.5 m. The bed material is mostly sandy and muddy. The shore’s slope is relatively steep, with an interval between 50° and 75°. At low water levels, the shoreline is typically 2–4 m above the water level (based on the Second Hungarian River Basin Management Plan).

Methods

There are several methods, approaches and models available for estimating hydrologic ecosystem services (HES) at landscape scale (Vigerstol and Aukema 2011). Hydrologic modelling tools are often specified along spatial scales (local or river basin scales) and in terms of temporality (quasi-permanent, event and process-based models). These methods differ in their data needs, in the level of user effort and in how sophistically different processes are described. The review published by EPA (US Environmental Protecting Agency 2018) clearly separates available tools into four classes (simple, mid-range, detailed/complex and modelling systems) with their classification following the complexity of the methods. Strengthening the descriptive structure of the models leads to a double effect: the uncertainty of the method used to describe the processes is reduced, but due to the increasing demand for data, the reliability of the data set is reduced (Chapra 1997).

The mainly GIS or cell-based ES models have common basics, and they are static models, as opposed to the tools described above, which are typically dynamic over time (Burkhard and Crossman 2013). Their purpose is to present typical ecosystem service values that occur during steady-state conditions of a longer period (from 10 to more year interval). There are also many differences, like required effort, integration and generalizability (Bagstad et al. 2013; Ochoa and Urbina-Cardona 2017).

We applied the InVEST model to map and quantify HES on the study area (Sharp et al. 2018). Three out of its modules were used to quantify and map the ES flood control, nutrient retention and erosion control in the Zala River Basin. The advantage of InVEST is that its data requirement is moderate and does not require in-depth expert knowledge. In principle, the results can be calibrated using global and site-specific local parameters and can be adapted to real-world data. In addition, the algorithm makes a number of significant simplifying assumptions in order to approximate environmental processes.

InVEST model has been compared in several studies to measured quantities and results of SWAT model, which is widely accepted and applied from a hydrologic point of view (Dennedy-Frank et al. 2016; Lüke and Hack 2017; Cong et al. 2020). The cell-based ES mapping tool estimated well the spatial distribution, and hot spot location of SWAT results in the above-cited studies.

In order to assess the capacity of the ecosystem to provide four selected HES, these ES were described by the difference between two model variants: outcomes of the vegetated land as it is at present (Actual-LULC) were compared with results of an un-vegetated (Bare Soil-LULC) model version. These two model set-ups differed only in their LULC conditions, and all other input files were unchanged. Influence of soil, topography, spatial distribution of precipitation and land use-specific parameters was included into modelling in order to reflect exposure to run-off/erosion processes. Flood control was assessed using the InVEST Seasonal Water Yield (SWY) module’s quickflow and baseflow results, while erosion control was assessed using the Sediment Delivery Ratio (SDR) tool with its sediment export result. For nutrient retention, both the export of total nitrogen (TN) and total phosphorus (TP) were taken into account, because the two substances are influenced by different processes (Chapra 1997), and their distribution will show a distinctly different picture (Table 1). For evaluating the effect of a potential adaptive land use change on the examined ES, we elaborated a future scenario based on aggregated scoring, which we describe step by step in the following sections. All simulations were based on the three-level hierarchical LULC map showing the actual vegetation, and the parameterization was also done for each of the 24 individual LULC classes (Table 3). All ES were evaluated with a 25 m spatial resolution. During the calculations, all our input data were fitted to the digital elevation model (DEM). The rasters with coarser (meteorological and soil hydrologic raster) and better (LULC raster) resolution were generated by resampling using the nearest neighbour method, while the same position of the cells was solved by masking with the DEM raster.

Table 1 Interpretation of the examined HES
Full size table

Flood control

Out of the 19 modules available in the InVEST model, the spatial variability of flood control service is best described by the SWY module. Using distributed parameter approach, the SWY estimates how much each terrain cell contributes to the total amount of the basin’s outflow through surface run-off (called quickflow) and baseflow. The model is suitable for describing longer periods (5, 10, 20 and more years of intervals), but it also provides an opportunity to show the variability between months and months under the chosen interval. The method is based on cell-level water balance and topography-based accumulation calculations. For a given cell, the Curve Number quantifies the proportion and infiltration and excess precipitation (NRCS 2004). Based on this, the model determines the amount of the surface run-off and reduces the precipitation by the actual evapotranspiration (ET) and quickflow, and local recharge is calculated. Summing up the local recharge along the topography-based flow paths leads to the baseflow values. The model uses the Curve Number (CN) and the crop coefficient (Kc) as land use–land cover-dependent parameters.

The SWY module was used to estimate the flood control ES in the Zala catchment for the 1981–2010 period. We interpreted the ES as the average yearly amount of water retained in each cell by the vegetation cover. In principle, it is possible to calibrate the model as well: the average discharge in the river basin over the period under study can be matched by the SWY estimated run-off by summing the quickflow and baseflow volumes upstream of a given cross section. The model can be fitted using land use- and soil-dependent input data and parameters. While we used recommendations (Allen et al. 1998; NRCS 2004) to set the CN and Kc parameters (Table 3), in case of the Actual-LULC variant we also carried out manual calibration of crop coefficients to achieve better model fitting. After that, the calculated discharge results were compared with measurements of annual mean discharge values and previous computation about the ET. The long-term average annual discharge of the whole pilot area and within that the sub-basins were adjusted to the measured discharges (based on the Second Hungarian River Basin Management Plan), and the results of the parameter adjustment were acceptable. The difference between the modelled and measured discharge values and the extent of the differences were quiet similar to the results of our previous study (Decsi and Kozma 2019). In addition, the spatial distribution of ET from the model results was compared with a previous analysis’ outcome carried out for Hungary (Szilagyi and Kovacs 2011).

Nutrient retention

The aim of the InVEST Nutrient Delivery Ratio (NDR) module is to map and quantify the non-point source emission of dissolved and sediment-bound nutrients (which are mainly total nitrogen and total phosphorus) and their transport to recipient water bodies. In general, nutrient emissions are channelled to recipient water bodies through different transport routes over the landscape (with surface run-off and subsurface leaking) as a result of rainfall. A portion of the nutrients are retained by the vegetation and the soil during run-off, which is influenced by topography and other processes. The NDR considers only the effects of terrain morphology and vegetation. Nutrient loads from the river basin are estimated in two steps: (i) at the cell-level with a simplified mass balance approach and then (ii) along the flow path using delivery ratios (Sharp et al. 2018) and connectivity index, which depends on hydrologic linkage between sources and sinks (Borselli et al. 2008). The nutrient retention of vegetation is described by a removal factor, which potentially reduces the specific emission from a given LULC category. This module is also static, and it is possible to characterize longer periods of use.

Emissions and their bound and dissolved ratios can be defined by LULC classes. For emissions, except for croplands and urban areas and pastures, values for atmospheric deposition were provided. In the case of croplands, we did the load estimation on county-level nitrogen budget data, based on the diffuse nutrient emission estimation methodology of the Second Hungarian River Basin Management Plan (Jolánkai et al. 2015). On paved surfaces, we defined twice the area-specific atmospheric deposition value, while in pasture areas it was defined as 1.5 times the nitrogen deposition, due to land use-specific excess emissions. In the case of phosphorus loads, we relied on the sample data sets of InVEST. The land use-dependent parameters of retention were used as recommended in the model documentation (Table 3).

One of the simplifications of the NDR tool is that it neglects the stages of nitrogen cycle and it considers watercourses as recipients, so it does not calculate with any processes in the water bodies. As long-term measured in-stream nutrient concentrations are available at only six sections, the estimation of diffuse nutrient loadings from the measured data is burdened with great uncertainty. Therefore, it is not possible to reliably calibrate the model results. The NDR was used to compute the total amount of retained and exported total nitrogen (TN) and total phosphorus (TP) at or from the watershed.

Erosion control

Using InVEST Sediment Delivery Ratio (SDR) module, it is possible to quantify erosion control ecosystem service (the amount of retained sediment by vegetation compared to non-vegetated conditions). The method can be used to estimate overland erosion, to track sedimentation and identify potential deposition and accumulation sites as well as quantities of retained sediment amount by vegetation.

The module is based on the revised universal soil loss equation (RUSLE1), but also takes into account the terrain morphology in addition for calculation of the potential soil loss (Benavidez et al. 2018). The annual typical soil loss per cell is determined using RUSLE1, and then, the soil loss yield is reduced by the Sediment Delivery Ratio based on the terrain, run-off characteristics of the unit (cell position along the flow path). The SDR index is a function of the connectivity index (IC), which is based on the position of the given cell and its location between the recipients and sources, as well as the types of downstream and upstream vegetation (Borselli et al. 2008).

As a result of the model, we get quantities of soil loss per cell (Sharp et al. 2018). The model considers units that can be assumed to be watercourses (their flow accumulation values are greater than 1000), and it does not calculate with sediment processes in the streams, so calibration is not possible directly.

To calculate the RUSLE1 equation at the cell level, we defined C (crop management) and P (support practice) factors per land use class, while R (rainfall erosivity) and K (soil erodibility factor) had to be defined in raster format. The rainfall erosivity factor was calculated using the method determined by Deumlich et al. (2006), based on annual rainfall values. Values for C and P factors were defined by using literature suggestions and the InVEST-sample database (Benavidez et al. 2018). Soil erodibility factors were determined by the textural classes from the DoSoReMI.hu database (Pásztor et al. 2017), based on the United States Department of Agriculture (USDA) texture classes (Table 3).

Aggregation

After the quantification of the four ES, we produced an aggregate result to represent the services under consideration. We applied a scoring system, likewise previously developed ones (Qiu and Turner 2013; Schulp et al. 2014; Willemen et al. 2018). We ranked each service on a scale of one to five (or bad to excellent) at the cell level. The classification system was as follows: for services with results of negative values (TP, TN retention and flood control), we separated the results of negative (less favourable results with the Actual-LULC than the Bare Soil-LULC model version) and positive values. The negative values (disservice areas with the Actual-LULC) were separated by its expected value, dividing the negative values into ‘bad’ and ‘weak’ scores. The positive values were divided into three parts along the 33rd and 66th percentiles, giving the ‘moderate’, ‘good’ and ‘excellent’ scores. In the case of erosion control, the class boundaries were defined by the 20th, 40th, 60th, 80th percentiles. The scoring system is summarized in Table 2. The final aggregated result—the value of HES delivery of each cell—was derived as the geometric mean of all four scores. Using this, low scores had a greater impact on the aggregated result than higher scores.

Table 2 Rules of the classification for a single ES (E: expected value, P: percentile of subscript, x: actual cell value of given ES, X: whole data set of given ES)
Full size table

Future land use scenarios

Land use planning is an interdisciplinary process influenced by number of spatio-temporal direct and indirect, as well as sample area-specific and economic factors (van Vliet et al. 2015; Harrison 2017). These drivers mostly depend on the disciplinary of the research (van Vliet et al. 2015; Plieninger et al. 2016). Nowadays, the most common method for this is land use (change) modelling, and there are several models with different approaches and disciplinary background (Verburg et al. 2004). Land use change studies are the most effective support for decision-makers and stakeholders when they cover the widest possible range of scenarios (Vallet et al. 2018).

The aim at this point of our study was to investigate the usability of the results of the InVEST model hydrologic ecosystem services in order to establish land use planning. In order to assess the impact of future land use change, we hypothesized a change directly based on this analysis. Based on the outcome of this evaluation, a theoretical land use change was applied for areas that appeared to be critical (which we defined as aggregate scores under 2). During the land use change, we built also on the simulation results of the current land cover conditions. We examined which land use category had reached the best aggregated rating, and then, the worst performing areas were replaced by the highest scoring land use type. After that, we classified the new results with the rules based on the original land use conditions’ aggregation rules. This way we could compare the results of two variants per ES and also the aggregated scores. The ratio of the hypothesized land use change was calculated per sub-basin in order to reflect on the distribution of ‘lost’ and ‘gained’ areas per LULC classes. This approach significantly simplifies the range of drives considered as well as the land use change strategy; also, it has technical limits (e.g. standalone cells, not feasible areas for afforestation) (Table 3).

Table 3 Summary of used databases
Full size table

Results

Flood control

As part of the parameter adjustment process, the actual evapotranspiration results of SWY were compared with the estimates of Szilagyi and Kovacs (2011). In 70% of the area, the difference in ET was less than 25% compared to previous calculations. In the remaining areas, SWY typically overestimated previously calculated ET. This may be because the period we examined covered a longer period, whereas the analysis of Szilagyi and Kovacs (2011) examined a shorter dry period. In the case of discharges, it was also difficult that the available data set’s sampling temporal interval overlapped only partially (20 years out of 30) with the period of our study. The available discharge values typically represent a more humid period. The modelled values showed a strong linear relationship with the available observed data (R2 = 0.91) on the 22 sub-basins. This resulted in a 15% lower outflow for the entire river basin. The tendency in residuals of ET and SWY is in line: the overestimation of ET explains the estimated lower value of discharge. In the light of the detailed uncertainties of the observed data sets, we accepted these results for use in this ES evaluation.

Based on the map presenting the results of flood control ES (Fig. 1a), the spatial distribution of meteorological conditions (precipitation, temperature and evapotranspiration) is generally reflected by the amount of water retained. So, as precipitation decreases from west to east, there is a decreasing tendency parallelly in flood control. Patches of water bodies (swamps, wetlands) and urban areas have a neutral effect, which can be explained by the fact that in the case of the Bare Soil-LULC variant these areas were left unchanged. As the processes of these land covers are typically more complex and influenced by incidental artificial factors, we decided to exclude these areas from flood control ES valuation. Environmental conditions like greater annual precipitation, more favourable LULC types (in terms of water management) and better soil hydrologic properties provide the best relative/specific water retention in the eastern part of the river basin: in some regions almost half of the precipitation was retained.

From the point of view of flood control, forests proved to be the most favourable (average value around 310 mm ha−1 year−1). Grasslands also performed well (almost 295 mm ha−1 year−1 on average), while the croplands showed a bit lower retention values (around 260 mm ha−1 year−1 on average). Impervious areas performed the worst. As it was expected, within their own land use category, green areas in urban environments (almost 290 mm ha−1 year−1) and among croplands the orchards and the complex cultivation patterns (almost 300 mm ha−1 year−1) have produced particularly good results.

Erosion control

The pattern of the erosion control service is moderately reflected in the terrain. The linear regression between the slope inclination and the amount of sediment retained has a moderate result (R2 = 0.32). Erosion control was mainly high in areas where the slope is relatively high (above 15%), and the variation in retention can be explained (i) with the position of the given cell along the flow path and (ii) the differences in land use (Fig. 2).

Fig. 2
figure 2

Spatial distribution of flood control (a) and erosion control (b) HES on the watershed

Full size image

The results suggest that land cover conditions are also dominant. Considering erosion control, forests turned out as the most efficient land use class (about 7 t ha−1 year−1 retention on average). The grasslands, like in case of flood control, performed well, with about 3.5 t ha−1 of sediment retained annually. Retention at urban, croplands and water bodies areas was estimated to be weaker. As it was presented above, within their own land use category, urban green areas (almost 4 t ha−1 year−1), orchards and the complex cultivation category (around 7 t ha−1 year−1) have performed better than their class peers.

Nutrient retention

The two nutrient components exhibit significantly different behaviours: although some areas may have negative retention in case of phosphorus (excess nutrient release compared to Bare Soil-LULC model version, purple patterns in Fig. 3b), this is more relevant for nitrogen considering either spatial extent or the measure of negative retention (purple patterns in Fig. 3a). About a quarter of the river basin area leads to nutrient emissions. This phenomenon is explained by the assumption made in land use-dependent input data: we hypothesized that croplands in the Bare Soil-LULC model variant do not receive external nutrient input, but in the case of the Actual-LULC significant organic or chemical fertilizer application can be expected for cultivated crops. Overall, croplands (primarily arable land) act as nutrient release sites.

Fig. 3
figure 3

Spatial distribution of nutrient retention ES (a TN retention, b TP retention)

Full size image

The different propagation characteristics are illustrated in Fig. 3. The effect described above can be properly observed: in the case of nitrogen, negative values (cells characterized by the resulting nutrient emission) occur in a significantly larger number and range, while for phosphorus nutrient retention (positive values) is clearly dominant.

Overall HES results by LULC categories

Figure 4 shows the main statistical indicators of the four ES calculated with respect to the first-level LULC categories. In the case of flood control ES, the grasslands and forests produced the best values as described above, but they also appeared to be sufficiently stable (with small interquartile ranges), as did the arable land, but somewhat weaker (Fig. 4a). Urban areas and water bodies led to the same median, but they appear to be more diverse. The reason for this is that within the third-level LULC categories there are significant differences between their collective categories (for urban: densely built-up vs. urban green surfaces with or without trees; for water bodies: reservoirs vs. marshes, swamps). In the case of flood control and TN retention, the distributions within the categories are generally balanced; on the other hand, TN retention of croplands and flood control of urban areas and water bodies had a wide range. Erosion control of forests has shown outstanding efficiency (Fig. 4b), while in case of nitrogen retention, as described above, croplands act as quasi-emitter/polluter in overall (Fig. 4c). During the mitigation of the other examined nutrient, the difference in land use classes was narrower, but the croplands were the weakest performing here as well (Fig. 4d).

Fig. 4
figure 4

Statistical distribution of the four ES (a flood control, b erosion control, c TN retention and d TP retention) with respect to first-level LULC classes (1: urban, 2: croplands, 3: grasslands, 4: forests, 5: water bodies)

Full size image

Aggregated

According to the results of the classification method based on the four services, 44% of the catchment were classified as ‘bad’ or ‘weak’. More than half of the basin area turned out ‘moderate’ or ‘good’, and only 4% performed ‘excellent’ (Fig. 5a). The areas with ‘bad’ and ‘weak’ scores were dominated by croplands. On the other hand, most of the land with ‘moderate’ or better scores had high ratio of forested areas and low rate of arable lands (Fig. 5a). Interestingly, urban areas also appeared more frequently in zones with ‘moderate’ or better grade. Grasslands (2–9%) and water bodies (1–3%) had a tight interval of coverage ratio within score categories (Table 4).

Fig. 5
figure 5

Hot spots based on ES aggregated evaluation with the actual (a) and the alternative land use–land cover conditions (b)

Full size image
Table 4 Distribution of the land use categories’ spatial extent within the score classes according to the present state (P) and with a hypothesized future (F) land use change, and differences of covered area by score categories
Full size table

Land use change

The change of land use affected slightly more than 80 km2, which is approximately 6% of the whole river basin. The average size of these patterns was 2.8 hectares, while the largest contiguous area covered 2200 hectares, and almost all the ‘bad’ classified areas were croplands (92%). These zones (their centroids) were within 500 m from the nearest watercourse in the 75% of the cases, while half of them were within 300 m to the closest stream.

As the forests received the best aggregated assessment, we hypothesized afforestation in places with ‘bad’ scores (aggregated values under 2). Our assumed land use change has brought a positive change for the Zala River Basin. One consequence of this is a shift in score distribution towards higher scores (Fig. 5b). The coverage of areas with ‘bad’ scores has almost dropped by 70% (Table 4). The number of ‘weak’ cells increased, and units with ‘moderate’ score also rose after the hypothesized land use change. Due to the alternative surface cover conditions, croplands’ ratings improved significantly, before the change they had owned an average aggregate score of 2.05, while with the alternative version it rose to almost 2.30 points (on the 1 to 5 scale). This effected the evaluation of forests, their aggregate score dropped with around 6% in average (Table 4). The increase in amount of nutrient retention was 75% for TN and 8% for TP on the whole catchment. For erosion control, the intervention had a neutral effect, while the flood control service also increased by almost 10%.

Figure 6 shows the ratio of afforestation per sub-basin and its effect on the regulating HES. Land use changes in roughly 6% of catchment’s area resulted in significant improvement of regulating ES primarily in nutrient retention. The average score of total nitrogen retention in these sub-basins increased by 8–25% in score (corresponding to 4–8 kg ha−1 year−1), while in the case of phosphorus it ranged between 5 and 10% increase (0.1–0.2 kg ha−1 year−1). Also for flood control, modelled LULC change resulted in an average 20 mm ha−1 year−1 improvement; however, its effect was smaller and less predictable (in some cases, a small land use change ratio resulted in a great increase in flood control ES, while in some the impact of greater changes was below expectations).

Fig. 6
figure 6

Effects of the hypothesized land use change on the ES per sub-basin (maximum, minimum: the highest or lowest change in average score out of the four examined ES on the whole watershed)

Full size image

Discussion

In overall, our methodology proved to be able to support decision-makers and stakeholders, as we have presented the spatial variability of several processes (both individually and aggregated) that basically influence river basin’s water quality and water quality trends and localized critical locations and also proposed an alternative to mitigate the negative impact of these areas. We tested a method to evaluate ES in an aggregated way, based on which it is possible to determine which land use performs best under which conditions. According to the cumulative assessment, the forests and grasslands performed best in the provisioning of HES (they were particularly effective compared to croplands). This efficiency was most evident for filtration (nutrient retention) HES that are mostly provided by forests, grasslands, green areas in the urban environment and wetlands.

Interestingly, urban areas produced better results than expected for all four ES we examined. The reason for this is that, according to the land use map (Tanács et al. 2019), green areas in nonrural environment are present in a high proportion (73% within their class). The reason for the marked difference between TN and TP is that the former leaks through the topsoil more easily, while the dominant fraction of TP tends to adsorb to particles in the surface covering layer, which leaves the area mainly through soil erosion. Overall, our results showed that a significant proportion of arable land can be considered a source of phosphorus and nitrogen emission.

In our analysis, implemented a simple aggregating system, using this we were able to identify the extreme areas (both excellent and bad) and finally to discover the primary cause of the poor performance. It was noticeable that those areas were typically classified as ‘bad’, which were under agricultural cultivation, and were typically located near watercourses (50% of ‘bad’ scored areas were within 200 m distance from stream network), and had less favourable terrain (around steep slopes) and soil conditions (grainy texture). As forests achieved the best valuation, we hypothesized land use change in these ‘bad’ areas by afforestation to consolidate or improve bad hydrologic services. Our hypothesized land use change was based on that planting trees is an increasingly popular civilian or government activity in the pursuit to reduce nutrient loads and improve water management and also mitigate carbon dioxide outputs and thus climate change. However, for actual implementation, modelling is not sufficient, and local guidelines for planting criteria need to be taken into account (e.g. Bozzay et al. 2019). In our hypothesized future land use model, slightly more than 6% of the river basin area had been changed that resulted in significant improvements in more than 70% of ‘bad’ classified areas (increase in aggregate score). In unfavourable cases where relief was present along with coarse-textured soils, forests have not solved the problem of erosion either. In these cases, it is worth extending the intervention to upstream and downstream areas, as well as creating a transition on the fringes of suggested afforestation and actual land use.

Trade-offs in ES occur when the provisioning of one service is increased, while the provisioning of another service decreases thereby (Raudsepp-Hearne et al. 2010). In our study, we did not assess any trade-offs directly, as we focused on modelling of four HES which usually do not influence each other in a negative way (Chapra 1997; Sutherland et al. 2018). While it was mostly cropland pixels that were rated in the aggregated HES worst, and thus selected to be turned to forest, the hypothesized changes increased forest-related ES (e.g. timber production, recreation; apart from HES) and decreased arable land-based ES (food production). The provisioning ES presumably decreased on the hypothesized afforested croplands, especially grain yields. Because in this study only regulating, but no provisioning ES were examined, it could be assumed that the crop yield would be decreased proportionally to the area of afforested land. According to that, we could expect a loss of about 15% grain yield and about 15% long-term increment of provisioning as well as other regulating ES of forests not examined here (e.g. carbon sequestration).

Due to the model simplifications by the processes’ description, the quantified results should not be interpreted primarily in absolute terms (e.g. exact amount of surface run-off, sediment and nutrient flows), but rather as an approximate estimation of spatial patterns, and those are typical of the area under study. On the other hand, InVEST performed reliable results in case of flood control, so in a simplified way, it considers the physical–biological background of the processes described. It is characterized by medium data and labour requirements and does not require any special in-depth expertise. At the same time, both the compilation of input map data and the choice of model parameters allow the utilization of expert knowledge.

Aggregating values from different ES is one way to deal with multiple aspects, which otherwise cannot be simply added up. Monetary valuation is another possibility that can help to unify a diversity of ‘goods’ in the widest sense, but this also has numerous difficulties and shortcomings (Bateman et al. 2011; Spangenberg and Settele 2016). Quantifying the importance of the single ES or ES bundles remains a task, where human weighting (according to the perceived importance of ES) is needed from stakeholders and/or from experts (as implemented, for example, by Harmáčková & Vačkář (2018)). Both scoring and aggregating are to a certain degree based on arbitrary decisions; however, the usefulness and scientific reproducibility of this process can be guaranteed by meticulous documentation (Czúcz et al. 2018; Schröter et al. 2018). In this study, we did not weight the individual ES in the final evaluation, but only aimed to assess where the problem is most pressing.

Considering the entire catchment area, a minor hypothesized afforestation resulted in large increase in the values of nutrient retention, especially in areas close to watercourses, while modelled land use changes on very steep hills or sandy soils produced sub-optimal results in flood and erosion control. The results indicate that interventions of relatively small spatial extent can produce significant results, when they occur at the site of critical areas. Knowledge of these critical focal points is essential in spatial planning, as the assessed processes influence the water quality of the confluent streams in general and of Lake Balaton in particular. Furthermore, when planning green infrastructure, these focal sites can be of great value for other spatial delineation tasks related to, for example, green infrastructure development or nature conservation measures.

Conclusion for future biology

Land use conflicts occurring at different scales, from regional to global, could be potentially resolved more effectively by objective analysis and modelling using an ecosystem services framework approach (Bagstad et al. 2013). In order to ensure the right quantity and quality of waters, modelling hydrologic ecosystem services should be an essential element for future decision-making. Several tools are available for this purpose, of which we tested a commonly applied method, the InVEST model for the first time in Hungary. Such relatively simple and easy-to-handle, spatially explicit models enable us to delineate areas where future change in land use might be potentially beneficial. Taking even more ecosystem services into account gives the opportunity to balance between more aspects, for example crop yield. Thus, present day decisions supported by our modelling results can influence the provisioning of future ecosystem services in a positive way.